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Abstract 

A roughly constant temperature over a wide range of densities is main- 
tained in molecular clouds through radiative heating and cooling. An 
isothermal equation of state is therefore frequently employed in molecular 
cloud simulations. However, the dynamical processes in molecular clouds 
include shock waves, expansion waves, cooling induced collapse and baro- 
clinic vorticity, all incompatible with the assumption of a purely isother- 
mal flow. Here, we incorporate an energy equation including all the im- 
portant heating and cooling rates and a simple chemical network into sim- 
ulations of three-dimensional, hydrodynamic, decaying turbulence. This 
allows us to test the accuracy of the isothermal assumption by directly 
comparing a model run with the modified energy equation to an isother- 
mal model. We compute an extreme case in which the initial turbulence 
is sufficiently strong to dissociate much of the gas and alter the specific 
heat ratio. The molecules then reform as the turbulence weakens. We 
track the true specific heat ratio as well as its effective value. We analyse 
power spectra, vorticity and shock structures, and discuss scaling laws for 
decaying turbulence. We derive some limitations to the isothermal ap- 
proximation for simulations of the interstellar medium using simple pro- 
jection techniques. Overall, even given the extreme conditions, we find 
that an isothermal flow provides an adequate physical and observational 
description of many properties. The main exceptions revealed here con- 
cern behaviour directly related to the high temperature zones behind the 
shock waves. 

*Email: gbp@phys.soton.ac.uk (GBP); m.d.smith@kent.ac.uk (MDS); mordecai@amnh.org 
(M-MML) 
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1 INTRODUCTION 



State-of-the-art numerical simulations of star forming regions take into account a 
small selection from a wide range of relevant physical processes. Of these, it has 
been found that supersonic turbulent motions combined with self-gravity and 
magnetic field reproduce o bservations of the mo rpholog y and dynamics of molec- 
ular clouds on all scales llTillev and Pudritd l2004t iBonnell et allEool l20n.lt 
IPadoan et afll2004 EooH lElmegreen and Scald see |2004 iMac Low and Klessenl 
see I2004L for reviews! . However, in order to reduce the computational con- 
straints, many results are based on simulations of an isothermal gas. This 
implies (1) highly radiative shocks, (2) suppression of the steepening of sound 
waves, (3) suppression of thermal instability and shock wave overstability, (4) 
no adiabatic cooling in expanding regions, (5) suppression of baroclinic vorticity 
and (6) helicity conservation. These assumptions are not always valid and could 
lead to erroneous results. In addition, one can examine neither the influence 
of the internal properties of the gas, such as its chemical composition, on the 
turbulent fluid flow nor the effect of the field of shock waves on the chemical 
composition. 

As one approach toward a more realistic scenario, we have performed numer- 
ical simulations of supersonic turbulence in dense, initially uniform, molecular 
materia l, with an equation of state determined by the chemical composition of 
the gas l)Pavlovski et all 120021 ; hereafter Paper I). We include a small but rele - 
vant chemical reaction network, based on the work of ISmith and Rosenl l|2003j) . 
We follow the time-dependent hydrogen chemistry, with C and O chemistry 
limited to reactions with ne utral atomic and molecu lar hydrogen, which gener- 
ated OH, CO and H2O fsee lSmith and Rosenl 120031 Appendix B). Our cooling 
function contains H2 ro- vibrational and dissociative cooling, CO and H2O ro- 
vibrational cooling, gas-grain, thermal bremsstrahlung and a steady-state ap- 
proximation to atomic cooling. That is, we follow shock-enhanced chemistry 
rather than the chemistry of the cold molecular gas. 

We demonstrated in Paper I, quite surprisingly, that even in the somewhat 
extreme case of high speed H2 dissociative turbulence, many integrated proper- 
ties of turbulence do not differ much from the corresponding properties in the 
isothermal regime. In particular, the rate of loss of energy in decaying turbulence 
obeys a similar law. Here, we examine non-isothermal turbulence in a broader 
context to better understand where to trust the isothermal approximation. 

Following the chemistry reveals properties not encountered in the isothermal 
simulations. Most importantly, we found that supersonic turbulence speeds up 
the refo rmation of hydroge n molecules by a factor of a few. This was devel- 
oped bv lSmith et alJ l|2004j) into a scenario for simultaneous rapid molecule and 
rapid molecular cloud formation out of atomic clouds. The accelerated chem- 
istry occurs in the dense layers behind shock waves where atomic collisional 
rates are enhanced. After a dynamical time, the reformed molecules are widely 
distributed as most of the gas has at some time passed through a shock. 

Here, we examine three dimensional decaying supersonic turbulence. We 
restrict our analysis to the high speed case of Paper I (see § |3 for description 
of the initial conditions; turbulence is chosen to be initially strong enough to 
destroy the molecules). The molecules subsequently reform and the field of 
shock waves then slowly decays within the molecular gas. Hereafter, we will 
use molecular turbulence (MT) as a shorthand for the run using the molecular 
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chemistry network data and isothermal turbulence (IT) for the run using the 
isothermal equation of state. 

In § 13.21 the power spectrum and cascade are analysed. Probability distri- 
bution functions (PDFs) for density, molecular density and velocity are pre- 
sented in § 0] The observational implications are considered in § Over 
the years, a number of different techniques have been pr oposed to in t erface 
analyses of numerical simulations with the observations llBrunt et all 120031: 
FIevcr and Schloera ll997tlOssenkopf and Mac Lowll2002tlMac' Low and Ossenkopl 
2000i IPadoan et aU Il998l) . To perform our comparison study, we apply here 
the same techniques to simulated observations derived from models based on 
both the chemical network described above and the isothermal assumption. 

2 The Simulations 

The gas is modelled within a three-dimensional box with periodic boundary 
conditions, simulated on a grid of 256 3 zones using a modified version of the gas 
dynamics code ZEUS-3D l|Stone and Normarl Il992j) . Paper I gives full details 
on the modifications to ZEUS-3D and on our initial set up. 

The initial state of the gas was chosen so that we could investigate a diverse 
range of conditions in the subsequent evolution. The gas is initially fully molec- 
ular. The imposed velocity field corresponds to that of the highest speed turbu- 
lence from Paper I, with a root-mean-square velocity of v rms = (60 kms -1 )t>6o- 
The turbulence is allowed to decay. In this manner, we can study the initial 
brief period within which shocks form and molecules dissociate. This is of order 
of the dynamical time 

tdyn = L/(kv lms ) = (15yr)Li 6 VQQ (k/S.B)' 1 , (1) 

for a box of size L = (10 16 cm)Li6 and a mean driving wavenumber k. This is 
followed by a period of molecule reformation of expected duration 

t Iei = 10 10 /(nT 1 ' 2 ) yr = (1060 yr)^ 1 ^ 72 , (2) 

for a hydrogen nucleon density of n = (10 6 cm~ 3 )ng and initial temperature 
T = (100 K)Xioo- This stage is considerably shortened by the turbulent com- 
pression. Finally, we follow an extended period of gradually decaying molecular 
turbulence. We set our numerical initial conditions to have Lie = n 6 = ^60 — 
7\oo = 1- 

It should be noted that the total cooli ng length behind a fast dissociative 
shock at the above density is ~ 10 cm IjSmith and Roserl |2p03) and is not 
resolved in the simulation, so we cannot predict quantities or instabilities in- 
herent to the shock physics and dynamics. The corresponding cooling time is 
approximately 

^cooi = (1 yr)ne 1 . (3) 

However, mass and momentum are still conserved, equilibrium chemistry is 
maintained, and the dissociation speed limits are approximately correct (Pa- 
per I). In the subsequent slower shocks within the molecular turbulence, the 
radiative layers are thicker and partly resolved. 

The temperature in IT was fixed at 100 K. All other parameters were un- 
changed from the MT case. That includes the initial kinetic energy spectrum 
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which is imposed arti ficially by perturbing the velocity field with a Gaussian 
random field following lMac Low et alJ l|l998j) . A narrow range of wavenumbers, 
3 < |k| < 4, is chosen for the velocity perturbation. 

It should be noted that the molecular simulation still excludes certain aspects 
of a molecular supe rsonic flow. The resolutio n precludes the appearance of the 
shock overstability l|Smith and Rosenl 1200^ . The zone size limits the degree 
of compression although high densities can occur when material accumulates 
within a zone. Both these effects can only be adequately modelled by decreasing 
the zone size by a factor of 100 or by introducing a strong magnetic field. 
Equilibrium chemistry, however, allows us to proceed and to follow some of the 
consequences of non-isothermal shock physics. 

3 POWER SPECTRA 

3.1 Definitions 

Power spectra of the fluid variables have long been used to characterise tur- 
bulence and so provide a standard means of comparing simulations. The total 
spectral power density is 

P{k) = ^ k )| 2 < ( 4 ) 

fc<|k|<fe+l 

where v(k) = T [v (r)] is the Fourier transform of the velocity field (and T is a 
Fourier operator). 

According to the Hclmholtz decomposition theorem, any differentiable vector 
field with divergence (V • v) and curl [V x v] can be split into compressional 
and solenoidal components. We thus write the velocity of the fluid as 

v (r) = v s (r) + v c (r) , (5) 

where 

[V x v c ] = 0, (V.v.) = 0. (6) 
The solenoidal and compressional power spectra are written as 

fc<|k|</c+l 

and 

P c (k)= J2 l*c(k)| 2 . (8) 

k<\k\<k+l 

The total power spectrum can also be expressed as P(k) = P s (k) + P c [k) since 
(v*.v c ) = 0. Note that we consider here angle averaged quantities, which is 
appropriate for isotropic turbulence (by isotropic we mean that the correlation 
length is smaller than the computational domain; see Paper I for details). 

Initially, the power spectrum of velocity also corresponds to the kinetic en- 
ergy distribution (or 'energy spectrum') since the density is initially uniform. 
The initial energy is divided one-third into compressive modes and two-thirds 
into solenoidal modes. 
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Figure 1: Evolution of the velocity power spectrum for IT case (thick grey lines) 
and MT (thin black lines) simulations. Data at three instances are displayed: 
upper lines: t — 30 yr; medium lines: t — 300 yr; lower lines: t — 600 yr. See 
also Fig. □ 

3.2 Total Power Spectra 

The spectral peak rapidly widens at the beginning of the simulations, as illus- 
trated by the curves in Fig. ^ for MT (thin black line) and IT (thick grey line). 
It develops into a wavenumber cascade with a maximum in the same initial 
range 3 < |k| < 4 for MT. The minimum is located at the high k end. 

At later times, the power is redistributed into a monotonic cascade with the 
maximum at k = 1 and minimum at k = 127 (see Figs.^l. For comparison, we 
indicate on the diagram the power-law Kolmogorov spectrum, A;" 11 / 3 , which is 
expected to occ ur even in decaying three dimensional incompressible turbulence 
l|Lesieurl Il997j) . In contrast, no clear power law develops in the spectra here 
as there is no continuing input of energy in the driving range in these decaying 
simulations and the power falls somewhat faster at high wavenumbers than at 
low wavenumbers in both cases (in this respect, the apparent 'convergence' of 
the curves at the high k end is illusory; see Table for the quantitative measure 
of the slope angle). In terms of the spectrum of shock waves, the increasing 
dominance of weak shocks preferentially dissipates the energy on small scales 
(see 8lO)l. 

The main result is that there is no significant difference between the evo- 
lutions of the total power spectra of MT and IT. Since the turbulence is not 
driven, the statistical properties of the fluid continue to change (as opposed to 
the dynamical equilibrium reached in driven turbulence simulations). 
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Figure 2: Evolution of the velocity power spectrum (ratio of compressional to 
solenoidal energy spectra) for MT on the left and IT on the right. Data at three 
instances are displayed: thick lines: t = 30 yr; thin lines: t — 300 yr; dashed 
lines: t = 600 yr. See also Fig. Q Typical amount of variability due to the 
limited number statistics can be assessed from the error bars plotted for the 
t = 600 yr IT case: values for large scales have smaller statistical significance 
as there are fewer large scale modes (the error bar for k = 1 isn't shown but is 
large since it is represented by a single value). 



Time MT IT 

Wyr' -2.91 ±0.06 -2.88 ±0.06 

300 yr -3.07 ±0.07 -3.07 ±0.07 

600 yr -3.16 ±0.07 -3.16 ±0.07 



Table 1: Coefficients, n, of linear fits to the total energy power spectra (logP oc 
n log k) for the two cases of molecular gas (MT) and gas with an isothermal 
equation of state (IT). The indices were measured over the interval k € [8,64] 
and demonstrate a systematic steepening of the curves with time. 
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3.3 Compressional and Solenoidal Power Spectra 

The power in compressional and solenoidal modes begins in a 1:2 ratio, as noted 
above. Intuitively, it might seem that the power ratio should scale with Mach 
number since a flow with zero Mach number is incompressible (all kinetic energy 
is in vortices) whi le a flow with a very large Mach number should be highly 
compressible ijElmegreen and Scald 12004). In our case, however, solenoidal and 
compressional modes are coupled and can exchange energy in both directions. 
Of special relevance here is how the evolution of vorticity, w = [V x v], is 
coupled with the divergence, (V • v), 

— + (v.V) w = (w.V) v - w (V.v) 
VI - x VP 



where p is the mass-density of the fluid, and P is pressure. The first term on the 
right hand side of equation © describes the vorticity generation due to tilting 
and stretching. The second term couples generation of the vorticity and diver- 
gence (if there is a negative divergence in the fluid, it will tend to shrink a fluid 
parcel, decreasing its moment of inertia, and, as a result, angular momentum 
conservation will give rise to vorticity), and the third term describes baroclinic 
vorticity generation. If turbulence is mo delled as barotropic or i sother mal, baro- 
clinic vorticity generation is suppressed (|Elmegreen and Scald . 120041) . 

As displayed in Figs. |3 the ratio of compressional to solenoidal energy does 
not remain constant across the range of scales, with large scales being more 
solenoidal and small scales slightly more compressive. This is not a surpris- 
ing result as compressive motions are associated with sho cks which are only a 
few z ones wide (as allowed by artificial viscosity (see, e.g.. lStone and Normarl 
[HH)) and thin shells from su bsequent r adiat i ve cooling. The shells break 
up due to various instabilities JyishniacL [l983; M ac Low and Normarl Il993t 
IVishniad.ll994HDgani et~ail Il996ft . " 

It is remarkable that small scale compressional motions develop much faster 
in the MT case than in the IT case (we thank M. Bate for poin ting this out). This 
effect can be partly due to additional thermal instabilities l|Smith and RosenL 
2003) that contribute to the rapid braking of shocks in the MT case, but is 
more likely to be associated with the higher sound speed in the dissociated gas 
that allows compressive waves to steepen faster. In the IT case, large scale 
shocks survive longer and compressional energy gradually shifts from large to 
small dissipative scales. Another possibility is additional baroclinic vorticity 
generation in the MT case, as follows from the analysis of the equation for 
vorticity evolution, Eq. 10. Statistical results for shocks and vorticity will be 
discussed in detail in § 14.41 

Once developed, the IT becomes more compressive at high wavenumbers 
than MT, remaining, however, predominantly solenoidal in both cases. The in- 
equality P c /P s < 1 has also been found in a simulation of driven superson ic mag- 
netised turbulence l|Boldvrev et all l2002HPadoan and Nordlundl l2002j) . Given 
that the total (solenoidal plus compressional) energy is the same in both simula- 
tions (see Fig. 1), we attribute this difference to baroclinic vorticity generation 
in MT. 
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A related result was also found in simulations of driven turbulenc e with 
a multiple power-law cooling function l|Vazauez-Semadeni et all ^996). The 
initial turbulence was purely solenoidal but was driven by purely compressible 
forcing. Vorticity production in this case was found to depend heavily on the 
additional terms (including Coriolis force on large scales, and Lorentz force on 
small scales). 

This interplay between modes was also noted bv lBalsara et al.l(|2004|) . There, 
even though the medium is driven by highly compressible motions, the kinetic 
energy is mainly concentrated in solenoidal rather than the compressible mo- 
tions. This behaviour was shown to arise from the interaction of strong shocks 
with each other and with the interstellar turbulence they self-consistently gen- 
erate. 



4 SIMULATED STRUCTURE 

4.1 Density PDFs 

The basic statistical property of a molecular cloud model or, indeed, observa- 
tional data is the distribution of density. While the observational picture is 
essentially two dimensional, our simulations provide detailed information about 
the density distribution in 3D. This distribution can be described in several 
ways. The fractional volume distribution is 

dv v = n v dp, n v = ^-^f, (io) 

V dp 

where dVy is the probability to find a value of density p£ [p,p + dp] and Vb is 
the total volume. Similarly, the fractional mass distribution is 

dv M = n M d P , o M = -L-^, (ii) 

M dp 

where M = (p)V~o is the total mass and (p) is the volume averaged density. 

To describe molecular turbulence, we introduce the fractional molecular 
mass, 

dP m = n m dp, n m = -±-^, (12) 

UJCo dp 

where 9Jt = (2/lA)fpdV is the mass of H2 (here, fn x 2m(H) is the mass of 
H2, and p = 1.4nm(H) is the density of the mixture: H, H2 and 10% of He) and 
97to is the total mass of H2 . 

The distributions given by equations (|10|) - (|12|l are normalised: 

n E = l, E = {V,M,£OT} (13) 

where S represents any index. Note that the Ox are not independent charac- 
teristics of the density distribution since it follows from equations (|10|l and (jl 1|> 
that 

M (p) = /rO v (p). (14) 
(P) 

Previous analyses of Oy and Om in simulations of compressible hydrody- 
namic supersonic turbulence have demonstrated that when the equation of 
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state is approximately isothermal, the density distributions are close to log- 
normal, i.e. that the logarithm of density has a probabil i ty distribution func- 
tion (PDF) that is a Gaussian llVazauez-Semadenl Il994t IPadoan et all Il997t 
Passot a,nd Yaznnez-Semade'nl Il998t IScalo et all ll99St lOstriker et all l200lh. 
Scalo et all l|l998h found that a log-normal density PDF should only occur at 



(1) low Mach numbers or (2) for any value of the Mach number when the poly- 
tropic index 7 is equal to unity (i.e. in the isothermal case). Otherwise, they 
found density PDFs to be well approximated by p ower-laws ov e r wid e regions 
of log(p) variation. The simulations performed by IScalo et all l|l998j) were of 
two-dimensional, driven, galactic turbulence, and included a variety of phys- 
ical processes expected on such scales, including self-gravity, magnetic fields, 
rotation, heating and cooling. 

It is not clear, however, whether the polytropic in dex (ratio of spe cific heats) 
alone is responsible for the shape of density PDFs. iLi et al.l (j2003) studied a 
number of cases of turbulence with an adiabatic equation of state with differ- 
ent polytropic indices. In the case of driven hydrodynamical turbulence, it was 
found that both mass and volume PDFs show imperfect log-normal distribu- 
tions. During dynamical evolution, the PDFs in the non-isothermal cases were 
found to develop only small deviations from Gaussian fits, while PDFs in the 
isothermal case (7 = 1) remain very well fitted by a Gaussi an. Simulations 
of thr ee-dimensional, supernova-driven, galactic turbulence bv iMac Low et all 
(2005) also show generally Gaussian results even th ough he a ting a nd cooling 
are included explicitly, consistent with the results of iKlessenl (|2000) who stud- 
ied the PDFs of isothermal, self-gravitating turbulence. However, he found that 
the moments of the density PDFs vary as collapse proceeds. Hence, self-gravity 
is also an important factor in shaping density PDFs. It is plausible that other 
physical processes can affect the shape of the density PDF. 

For parameterisation of PDFs with log-normal distributions, it is useful to 
define a dimensionless variable x — m (p/(/°))- Substituting into 

^=n E (p) (is) 

dp 

yields 

Vn(x)- = ^(P)- (16) 
P 

If Qs(p) is a log-normal distribution, then Vs{x) should be the normal (Gaus- 
sian) distribution 



Vt.{x) = -7= — exp 

V27T<7£ 



2o* 



(17) 



Following the work of lLi et all l)2003l) . we relate Vm(x) an d Vv{x)i using equa- 
tions HHJ and JI3||: 

Vm(x) = e x V v (x)- (18) 

For the functional dependence given by equation (|17(1 to be valid for both PDFs, 
parameters of the normal distributions have to satisfy the following relations due 
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0.0 0.1 0.2 0.3 0.4 0.5 
fraction, <f> 

Figure 3: The average specific heat ra- 
tio (7) employed by the code as a func- 
tion of the average molecular fraction 
(/) in a warm molecular/atomic gas 
mixture (solid line) ; the isothermal gas 
case 7 = 1 is also shown (dashed line). 



10 100 

time [yr] 

Figure 4: Long dashed line: (7) the 
average specific heat ratio as a func- 
tion of time in the MT simulation; 
Continuous line: 7 P the effective spe- 
cific heat ratio as a function of time; 
short dashed line: the isothermal value 
(1.0) 



to Eq. dEJ: 

cry = <r M = c, (19) 
XM=Xv + <J 2 , (20) 

XV = — 2", > XM = y , (21) 

as noted bv lOstriker et all l)200l|) . 

In terms of the PDFs, the distributions ilgji(p) an( i ^v(p) are related by 

^ot = j-J-rSlvi (22) 
ifp) 

where / is molecular fraction (fn is abundance of H2), (fp) = (1.4/2)9Jt/V. 
Using equations Ij22(l . Ijl4|l and Ijl6(l we can relate the normal distributions for 
mass and molecular mass through 

4.2 The Effective Specific Heat Ratio 

The actual ratio of specific heats, 7, as calculated for each zone in our MT 
simulation, is a monotonic function of the molecular fraction, /, 

7 =(5.5-3/)/(3.3-/), (24) 

(including 10% helium). As shown in Fig.0 7 may vary between 5/3 and 10/7. 

We can define a spatially-averaged specific heat ratio by substituting / — ► 
(/) in Eq. 124|) . With these definitions, (7) does not depend on the gas exci- 
tation. Hence, strictly speaking, it can serve as a good approximation only for 
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Figure 5: Scatter plots: logarithm of pressure vs logarithm of density. Grey 
points correspond to the time t = 10 yr, and black points t = 600 yr; linear 
fitting log(P) oc log(jd) has a coefficient close to 1.0 at later times, see Pig. SI 

a gas at a high temperature, to ensure that the rotational degrees of freedom 
are fully excited and the classical approximation for specific heats, which we 
employed in our model, remains valid. For molecular hydrogen, this is true for 
temperatures T > 85 K. As shown in Fig. 21 (7) varies considerably during the 
simulation. 

To compare to the isothermal simulation, we define a polytropic index j p . 
This is the effective specific heat ratio including radiative cooling. We define it 
as the linear fit coefficient for the dependency 

log(P) = 7p log(p) + const, (25) 

where P is pressure and p is total mass density. In Fig0 we display the P — p 
data at two times. At early times (grey cells), the dependency displays a great 
deal of scatter but is effectively bound between two j p — 1 lines which, overall, 
results in the linear fitted value being close to 2, see Fig.^] The upper line -f. p — 1 
is formed by the highest temperature regions; this temperature corresponds to 
~ 8000K, the post-shock temperature following rapid atomic cooling but before 
molecular cooling has become effective. The lower line is a reminder of the 
original T = 100 K isothermal conditions. At later times (black cells in Fig. 
the scatter is reduced, and it indicates an effective isothermality of the gas. 

4.3 PDFs: analysis 

In fully molecular gas f(p) = (/) = 0.5 and Vm = Pm- Hence, differences 
between these two distributions are only going to be noticeable when strong 
shocks dissociate a considerable amount of H2. As shown in Fig.0J this occurs in 
the first 100 yr of the simulation. Over the next 500 yr, owing to the reformation 
of molecules, we can expect large deviations of P<%i from Pm (see also Fig. 4 of 




11 



Paper I) . Finally, when dissociative shocks disappear and the molecular fraction 
distribution flattens, Vm and Vm should converge. 

The volume-weighted density PDFs in the MT run are found to be close to 
log-normal, as shown on the series of plots in Fig. (left panel). Volume- 
weighted PDFs of density for Burgers' turbulence (completely shock domi- 
nated turbulence) have been analytically s hown to have power-law asymptotics 
l|Frisch and BeckffeOOlHFrisch et aill200lj) . It is clear from the error bars on the 
plots that we can't make a definitive statement about the asymptotic behaviour 
due to the lack of numerical resolution, but our results do not show clear evi- 
dence for this behavior. Numerical observations of PDF asym ptotics require ex- 
tremely hig h resolution even in one-dimensional cases (see, e.g.. lGotoh and Kraichnarl 
Il993l Il99q) . We note however, that, during the early hypersonic evolution of 
both the IT and MT cases, the mass- weighted PDFs display clear power-law 
tails at low densities. We remain uncertain how to interpret this result. 

To determine how much the deviation from a log-normal distribution is in- 
fluenced by the high and variable value of the specific heat ratio and the in- 
troduction of explicit heating and cooling, we contrast these results with PDFs 
from IT (right panel of Fig. |BJ , which has exactly the same initial physical con- 
ditions and numerical properties. The PDFs correspond to the same moments 
of time. Both data sets suggest that deviations from log-normal distributions 
are more significant at earlier times. As turbulence decays, the distributions 
become more consistent with log-normal properties (see equations (|TT?)l — l|2"T|l 
). We find that the PDF of density from the MT simulation is actually more 
consistent with a log- normal distribution than that from IT. 

The evolution of the first four moments (mean variance of., skewness, 
/fe, kurtosis, <?s) of the PDFs are displayed in Fig.0 The variance, <j v , exhibits 
a large deviation from <j m early on during both simulations. However, in the MT 
case the deviation is less significant. In both cases, a v rises sharply when shocks 
form and distort the initially homogeneous density field. At the beginning of the 
simulation, in molecular gas the density contrast is less than in isothermal gas 
because 7 > 1. As shown in Fig0] a significant hot gas phase is produced. This 
can qualitatively explain smaller values of Oy in a molecular turbulent flow. 
The variances a\j and are approximately constant during the simulations. 
In the MT simulation, parameters of mass and molecular mass distributions 
converge, as expected from the above analysis. 

To measure deviation of the log-distribution from gaussianity we also com- 
pute higher order moments — skewness, 0, and kurtosis, <?, 

? = 9=~o, (26) 

where 

»n = ((x-(x)) n ) (27) 

is a central moment of n-th order. For a normal distribution, skewness, a degree 
of asymmetry of the distribution, (3 = 0, and kurtosis, a degree of peakedness 
of the distribution, g = 3. For an exponential PDF the corresponding values 
are (3 = 2, g = 6. As seen from plots on Fig. El a significant deviation from 
Gaussian shape occurs only at earlier times. However, skewness remains nega- 
tive throughout, indicating longer 'tail' on the side of low densities. It can be 



12 





Figure 6: Evolution of PDFs of the density distribution during the simulation 
of decaying, high-speed, molecular turbulence including heating and cooling 
(MT, left panel) and isothermal turbulence (IT, right panel). Note power- 
law wings of the mass-weighted PDFs, and more consistent log-normal shape 
of the volume- weighted PDFs. To demonstrate statistical significance of the 
data points we plot error bars for every 20th point of each PDF curve. When 
statistical significance of a PDF point is small (i.e., it is derived from a single 
data value) the curve is traced with a small horizontal dashes (see extreme of 
the wings of the mass- weighted PDFs). 
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Figure 7: First moment (mean) %s, second moment (variance) <7 S , third moment 
(skewness) /3s, and fourth moment (kurtosis) <?£ of Vs as a function of time 
in the MT (left) and IT (right) simulations. The volume mean, \v is plotted 
with inverse sign. A Gaussian distribution has skewness, (3 — 0.0, and kurtosis, 
g = 3, as indicated by the dotted lines. 

attributed to the fact that highest values of density are limited by the high- 
est Mach number, whereas the lowest density does not depend on properties of 
individual shocks directly. 

These data suggest that (7) only slightly influences the PDFs of the den- 
sity distributions. Although these results are more generally applicable than 
specifically to large-scale, cold molecular clouds, where we would expect values 
of polytropic index 7 w 1, it is remarkable that an entirely different equation of 
state leads to qualitatively and quantitatively very similar statistical properties. 

4.4 Velocity: vorticity and shock fields 

As remarked above, the assumption of isothermal or isentropic flow stops baro- 
clinic vorticity creation. It also forces helicity, to = ((v • w)), to be conserved in 
an unphysical way in general compressible flows. In this sub-section, we analyse 
in detail the differences in the amount of vorticity and shocks in the MT and 
IT cases. 

The vorticity distribution is displayed in Fig. In isotropic, incompress- 
ible, turbulent flow, high vorticity exists in thin coherent tubes, as found in di- 
rect numerical simulations (D NS) with high resolution (see, e.g.. lLesieur et all 
l2008t iTiton and Cadotl [2003ft . The vorticity tubes have a diameter of a few 
Kolmogorov dissipative scales and a length of the order of the integral scale 
(i.e. the driving scale). Centres of the vorticity tubes form one dimensional 
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Figure 8: Left cube: Surfaces of isovorticity, |[V x v| = const at 10 yr after 
the start of the MT simulation. Shocks have just formed and large structures 
are present. The black upper face of the cube contains a contour plot of the 
top slice of vorticity. Right cube: surfaces of isovorticity at the later time of 
300 yr when the shocks have fragmented. 

filament-like structures that are the centres of energy dissipation of the flow. 

In comp ressible flow, on th e other hand, vorticity structures form sheets 
and spirals l|Porter et all l2002h . In our simulations of decaying supersonic tur- 
bulence, we find that the vortex structures (isosurfaces of high vorticity) form 
"open" surfaces. The size of these surfaces is closely correlated with the size 
of shock surfaces (Fig. Isosurfaces of strong vorticity become smaller with 
time, with a characteristic size of order of the shock size, as displayed in Figs. [3] 
and El 

Statistical analysis of the vortex structures reveals that vorticity components 
as well as velocity gradients have PDFs of exponential type, i.e., proportional 
to exp{— \X\}, in incompressible turbulent flow. This is instead of the Gaussian 
distributions proportional to exp{— \X\ 2 } that would be expected for statisti- 
cally independent structures due to the Central Limit Theorem. In laboratory 
experime nts and DNS, regions of high vorticity correlate with regions of low 
pressure (|LesieurLll997l p. 196). The exponential form of the PDF of vorticity is 
a manifestation of the intermittency of energy dissipation. If dissipative struc- 
tures were distributed homogeneously, the PDF would be given by a Dirac delta 
function, or a narrow Gaussian distribution in a DNS. Therefore, a deviation of 
the PDF of dissipative structures from the Gaussian shape is a direct indication 
of the intermittency of turbulence. 

To study PDFs of shock waves we define the function of convergence, conv(w), 

f(V-v), when(V-v)<0, 
conv(v) = <\ 7 ' \ ' (28) 

I 0, otherwise. 

A simple statistical analysis of the converging regions and vorticity reveals the 
expected exponential PDFs for both fields fFig. H0 | ). Exponential law s for PDFs 
related to convergence were reported before bv lSmith et alJ (|2000bl) . where the 
distribution of shock jump velocities in a single direction were found to be 



15 





50 100 150 700 250 



Figure 9: Contour plots from the MT case of negative velocity divergence (black) 
and contours of the vorticity field (grey) . Note that the vorticity is concentrated 
around the shocks. The plots display random slices from the data cubes at 
t = 10 yr (left) and t = 300 yr (right). 



approximated by an exponential: 



dN = number of zones with Avi £ [Avi , Avi + e] 
oc exp{— CAvi}, 



(29) 



where C is a positive constant, and e is a small discretisation parameter. 

We find that the PDF of conv(v), which is essentially a 3D generalisation of 
Smith et al.'s one dimensional jump velocity, has the same exponential depen- 
dence, 

dN 

-=3 cx exp{-/3|conv(u)|}, (30) 

where dN is the number of zones with conv(w) €E [conv(t> ), conv(w) + e] and N 3 
is the total number of zones (see Fig. I10|) . We have chosen the discretisation 
parameter e to create 500 bins for our histograms. 

Vorticity PDFs are studied using the modulus of curl, |curl(i>)| = | [V x v]|. 
We find that the PDF of curl behaves very similarly to the PDF of convergence 
but with different parameters: 



^ oc exp{-a|curl(w)|}. 



(31) 



The coefficients a and f3 change as the turbulence decays. Both PDFs be- 
come progressively steeper with time, indicating that strong shocks and vorticity 
disappear. However, we find that the PDFs of convergence and vorticity do not 
change with time if we compensate for the average convergence and curl decay. 
The coefficients a and (3 of the distributions, 



diV , |curl(u)| , 

— —t cx exp^a-7i 77 ... k 

N 3 Fl (|curh»|) J 

dN conv(v) , 

— oc exp{/3- — }, 

iv J (conv(w)) 



(32) 
(33) 
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Figure 10: Left panel: histograms of |curl(ii)| at different times (MT simula- 
tion). Right panel: histograms of conv(w) at different times (MT simulation). 
The histograms correspond to the following simulation times: 10, 20, 30, 50, 
100, 150, 200, 300, 600 [yr], with the direction of the time sequence marked by 
the arrows on the figures i.e. the distributions steepen. 

MT IT 

a 2.63 ± 0.05 2.61 ± 0.06 
2.48 ±0.05 2.30 ±0.05 

Table 2: Coefficients of compensated PDFs of convergence (/3) and module of 
curl (a). The error estimates are derived from the statistical variance of the 
distributions of the coefficients derived from different moments in time). 



do not significantly change during the simulation. The values of the coefficients 
for MT and IT simulations are given in Table El These values may be significant 
constants characterising decaying turbulence. 

The data in Tabic [21 suggest that the compensated PDF of vorticity is not 
influenced strongly by a different equation of state, whereas the compensated 
PDF of convergence is steeper in the MT case. A shallower PDF of convergence 
implies the existence of larger quantities of strong shocks in the IT case, which is 
not surprising. Shock jump conditions imply that in the case when the ratio of 
specific heats 7 > 1, the velocity difference across zones of the shock is smaller 
than in the case when 7 = 1. This may qualitatively explain the steeper PDF 
of convergence in molecular turbulence. 

The PDF of the velocity itself is not very well understood l|Elmegreen and Scald 

l2004h . The usual argument is that the velocity PDF should be Gaussian is based 
on the central-limit theorem applied to Fourier components of the velocity field. 
This argument, however, does not take into account correlations which are im- 
portant for energy transfer among scales. In fact, the veloci ty PDF must p ossess 
non-zero skewness for the transfer to be possible (see, e.g.. (Lesieuit Il997h . 

Analytical studies of Burgers' (shock-dominated) turbulence foun d that the 
gradient of the velocity has power law asymptotics (same as density, see Frisch and I 
(200l|). We do not see this effect in our simulations, perhaps because of our 
limited numerical resolution (see also t !4.3l) . 

In our simulations, the initial distribution of the velocities has imposed Gaus- 
sian statistics. We find that deviations from the Gaussian distribution do occur 
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Figure 11: PDFs of the velocity magnitude during the MT simulation. The peak 
shifts toward smaller values as the turbulence decays. Plotted curves correspond 
to the same moments of time as curves on Figs. & I1UI The arrow points in 
the direction of time (velocities decrease). Dotted lines show Maxwellian fits. 



but they are minor. To illustrate this, in Fig.^Jwe plot the PDFs of the velocity 
magnitude, v = \/ x 2 + y 2 + z 2 , which takes the form of a Maxwell distribution, 
i.e. exp{— Av 2 }v 2 7 as expected. 



5 SIMULATED OBSERVATIONS 

5.1 Background 

Observational predictions from numerical results have generally in volved the 
post-processing of isothermal flows to generate s i mulated observations llFalg aroue 
1994HPadoan et allll998ULazarian et alll200lHOssenkopf and Mac Lowt T 



Brunt et all 2003^1 . We can validate such an approach by comparing synthetic 



maps from isothermal models to similar maps from the molecular models. We 
note that the compact and dense regime taken in our simulations is not directly 
comparable to existing observations of molecular clouds but rather corresponds 
to dense cores. 

Interstellar clouds appear to f ollow a scaling relation first inferred from an 
analysis of observational data bv iLarsonl l)l98l|i . Larson used 13 CO (optically 
thin line) maps of different molecular clouds to establish a size-linewidth relation 

a v (x\ b , 6 = 0.5 ±0.2, (34) 

where a v is the velocity dispersion and A is the size of the region. 

A theoretical explanation for this scaling relation can be inferred from the 
following argument. The velocity dispersion, a v , can be approximated by the 
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square root of the second order velocity structure function, 
S 2 = ((v(r) - v(r')) 2 ) oc A ?2 , A = |r - r'|. 



(35) 



The value of C2 is dete r mined from the model of driven supersonic turbulence 
proposed bv iBoldvrevI l|2002h . based on the dimensionality of the dissipative 
structures. Boldyrev's theory suggests the following scaling law: 

a v = A 37 , (36) 

which is in good agreement with the above cited observational value. Pure 
K41 scaling gives smaller values of the exponent (b « 0.33). However, in the 
case of decaying turbulence Boldyrev's scaling laws do not apply. As we have 
demonstrated in Section^ decaying turbulence results in a shallower power law 
energy spectrum than follows from K41 theory. To estimate the exponent b 
in this case we can use the connection between velocity and power spectrum 
statistics, 

Pi(A)ocfc- m , (37) 
S 2 oc\ m -\ (38) 

where P\{k) is the (one dim ensiona l ) pow er spectrum; using the notation of 
section [31 n = m + 2 (see, e.g. jFrischLfl995^ . Using values of n given in Tabled 
to « —0.88 ... — 1.16, which implies, b = —0.06 . . . 0.08. Hence, we can predict 
that a v should not exhibit strong scaling in the case of decaying turbulence. 

Deduction of Larson's scaling law from the observations is not a straightfor- 
ward task. The width of a spectral line, broadened by the turbulent motions, 
is essentially the width of a composite line, contaminated by all the gas mo- 
tions along the line of sight. How the width of such a composite line should 
depend on the size of the sampled region is not obvious. If we sample a v from 
the regions that appear to be cores (i.e. continuous prominent patches of ra- 
diation intensity), sorting the a v measurement by the size of such cores, we 
might stand a better chance of recovering true velocity scaling. This is because 
such cores may actually be coherent 3D structures, al though there is no way we 
can protect our statistics from unavoidable mistakes llVazauez-Semadeni et all 
Il997t IPichardo et all Eooot iBallesteros-Paredes and Mac Lowtl2002|) . We have 
no means of distinguishing between "false" cores — column density enhance- 
ments produced by overlapping, unrelated density enhancements assembled by 
chance along the line of sight — and true cores — compact 3D objects. In 
regions of active star formation where self-gravity effects are important, we can 
hope that the number of such mis-identifications will be effectively minimised. 

In the case we are studying, there is no boosting of contrast due to self- 
gravity but there are additional contrast parameters: the chemical composition 
and temperature distributions. Indeed, our simulations provide detailed pictures 
of the H2 and temperature distributions. Non- uniform chemical tracers may 
provide additional contrasts to the maps. What sort of regions can chemical 
tracers emphasise? As we found out in Paper I, the distribution of molecules in 
decaying turbulence is rather homogeneous, and there is no evidence that the 
H2 distribution alone will highlight 3D structures within the simulated region. 
However, the contrast in temperatures is far more distinct. Most of the gas in 
the simulated domain stays cool due to the efficient cooling; only sufficiently 
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Figure 12: Histogram of the temperature distribution at t — 60 yr. The his- 
togram has 1000 bins with size 6T = 16 K. A large fraction of the gas has 
temperature close to T = 100 K. The average temperature is w 630 K. Data 
taken from the MT simulation. 

strong shocks can temporarily heat up the gas in certain regions, up to several 
thousand degrees. Hence, if we select a tracer sensitive to the temperature, we 
would be able to construct maps of regions with different values of temperature. 

Temperature sensitive tracer maps may help answer another question rele- 
vant to the present study. Is detailed information about molecule and temper- 
ature distributions crucial for predicting emission, or is the zero-order approxi- 
mation of isothermality sufficient? 

5.2 Scaled projection of CO emission 

The gas in our simulation is relatively warm with an average temperature > 
100 K (Fig. 112(1 . Rotati onal emission lines of 12 CO are a suitable tracer in 
this temperature range l|McKee et all Il982). Emission from the transitions 
between states with different rotational quantum number J (only J — ► J — 1 
transitions are allowed due to the quantum selection rule) will highlight regions 
with different temperatures. The choice of this particular tracer is reasonable 
because our molecular code actually computes the equilibrium abundance of CO 
to account for the cooling through rotational and vibrational CO emi ssion. 

We take the non-LTE formulae presented bv iMcKee et al.l (^82). We as- 
sume that each emission line is optically thin and, therefore, the intensity of the 
emission (the "brightness" of the pixel on the emission map) is proportional to 
the column density of CO, 

emission at J level oc (nj (CO)) jAhj (39) 

where (nj(CO))i is the average number density of CO at level J in the cell of 
size (Ahi) 3 . 

The simplest way to construct a synthetic map of CO emission is to map 
the column density, as prescribed by equation (|39|) . Following the work of 
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lOstriker et al . (2001), we identify a clump as a region of contrast (ROC), which 
is regarded as such if its CO column density is at least a factor of c r larger 
than the mean column density of CO for the entire map. This procedure is 
somewhat similar to the procedure of subtracting "background" in the reduction 
of observational data. In our current study, we fix c r — 1, i.e. all regions with 
column density larger than the mean column density are identified as ROCs, 
and we constrain our analysis to the data we infer from the ROCs only. 

To deduce the relationship between line width and size from the synthetic 
maps, we need an algorithm for the detection of ROCs of different sizes. One 
means of achieving this is by constructing "blurred" maps of lower resolution 
by rescaling the original map of CO column density by a certain factor, as 
proposed bv lOstriker et al.1 l)200l|) . We rescale the original 256 2 grid by a fac- 
tor of 2 s , where s — 0, 1, 2, . . . , 7. This results in eight maps with resolutions 
2 2 , 4 2 , 8 2 , . . . , 256 2 , which we can employ for ROC identification. 

To construct CO emission maps from the IT simulation with T = 100 K 
we need to determine an appropriate abundance of CO (in each computational 
zone) without having information about the H2 abundance. The chemistry algo- 
rithm implemented in our version of ZEUS-3D implies that, for this temperature 
all free gas-phase carbon is in the form of CO for all plausible mass densities and 
H2 fractions. The total abundance of C is fixed in the code, f(C) = 2.0 x 10~ 4 . 
Hence, in the isothermal case we take the distribution of CO as homogeneous 
with the fixed value of 2.0 x 10~ 4 . 

We select the time when the largest fraction of the molecular gas is at T = 
100 K (see the temperature histogram in Fig. ll2f> . The actual shapes of lines and 
maxima in CO emission will depend on the underlying mass density distribution 
and velocity along the line of sight. 

To be able to compare the morphological features directly we use density 
and velocity data from the MT simulation in both cases. The actual distribu- 
tions of CO and temperature are employed for the MT maps, and homogeneous 
distributions of CO and temperature for the IT maps. We will refer to such 
isothermal data as quasi-isothermal. As we will show, the results of the analysis 
of the quasi-isothermal maps are very similar to the results of the analysis with 
the actual isothermal simulation data. This is not surprising giving that the 
statistical properties of the density and velocity fields in both simulations are 
very similar, as established in §§ 14.41 and 14. II 

For each ROC, we compute the mass- weighted dispersion of the line-of- 
sight velocity a v , which represents the line width for a region of projected area 
[(L/N)/2 S ] 2 . We record M(COj) within each ROC, and the virial parameter, 

So-l[{L,N)/n 
GM(COj) ' K ' 

where M(COj) is the mass of the COj (CO resp onsible for J — > J — 1 emis- 
sion) and G is Newton's gravitational constant l|Ostriker et "all. l200lh . The 
parameter a is directly proportional to the ratio of kinetic energy to gravita- 
tional energy, and represents a measure of the relative importance of gravity 
llMcKee and Zweibel 11993 lOstriker et all 1200 1|) . 
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5.3 Analysis of the emission maps 

We now present the results of the synthetic observations. For the data sets 
discussed in section 15.21 we have constructed CO (4-3) emission maps that 
highlight the gas at temperatures around 100 K (for gas with average density 
of 10 6 cm~ 3 ), and C O (20-19), which traces much hotter gas of about 1000 K 
i McKee et allll982|) . 



CO (4-3) emission maps for the MT data and the quasi-isothermal data 
are presented in Fig. EH The similarity of these maps is quite remarkable. 
Differences are difficult to detect. This supports the claim made in Paper I that 
isothermal simulations are suitable for mapping molecular emission. 

An analysis of the statistical properties of these maps also reveals no signifi- 
cant discrepancies. In Figs. ll4H17l we display scatter plots of the parameters of 
the ROCs which were identified on the MT maps at eight different resolutions. 
Each point or small dash on the plots corresponds to the parameters of a ROC. 
Although there is a general trend in mean dependences, there is plenty of scatter 
in all the plots. The solid lines represent least-square linear fits of the form 

dlop|A) = 
dlog A 

where A = {L/N)/2 S is the size of the ROC; 

dlog MA) =6 ,. (42) 

- c; (43) 

= d. (44) 

The parallel dashed lines in the figures mark the la deviations for the fits. The 
values iflTJ) - (|4*4*jl deduced from the MT data are as follows, 

b = (2.6 ± 0.4) x 10~ 2 , b' = (1.7 ± 0.2) x 10~ 2 , 

^ V ' (45) 

c = -(4.94 ± 0.04) x 10"\ d = 1.984 ± 0.003. 

The values b, b' , c and d deduced from the map constructed from the quasi- 
isothermal data are the same within the error estimates: 

b = (2.7 ± 0.4) x 10 -2 , b' = (1.7 ± 0.2) x 10~ 2 , 

; V ^ (46) 

c = -(4.86 ± 0.04) x 10" 1 , d = 1.985 ± 0.003. 

As noted in § 15.21 the statistical properties of the emission map constructed 
using density and velocity distributions from the actual IT simulations are vir- 
tually the same. The value of parameters in that case are 

b = (1.7 ± 0.4) x 10~ 2 , b' = (1.1 ± 0.2) x 10~ 2 , 
c = -(4.99 ± 0.04) x 10"\ d = 1.985 ± 0.002. 

Such a good agreement in the morphology and statistical properties of the 
maps is due to the very homogeneous distribution of molecules in the gas. At 
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Figure 13: Maps of 12 C0 J = 4 — * 3 emission constructed from the MT 
data (upper panel), and the same emission map constructed from the quasi- 
isothermal turbulence data (lower panel), as described in section 1^1 
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Figure 14: Sigma vs. size of ROC 
for the J = 4 — > 3 molecular tur- 
bulence map. For reference, c s is 
the mass-weighted average speed of 
sound. 
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Figure 16: Alpha parameter vs. 
mass of ROC for J = 4 — > 3 molec- 
ular turbulence map 
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Figure 15: Sigma vs. mass of ROC, 
for J = 4 — ► 3 molecular turbulence 
map. For reference, c s is the mass- 
weighted average speed of sound. 
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Figure 17: Mass of ROC vs. size 
of ROC, for J = 4 — » 3 molecular 
turbulence map 
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this stage of molecular turbulence evolution, t — 60 yr, there are numerous 
dissociative shocks, and the mean characteristics would suggest that molecular 
turbulence should be very different from isothermal turbulence in many respects 
(at t = 60 yr the average molecular hydrogen fraction is (/) =0.1, and average 
temperature is (T) = 630 K). Despite this, the analysis of actual molecular 
emission maps provides evidence to the contrary. 

Furthermore, the shapes of spectral lines in both cases are very similar. The 
line shapes are very roughly Gaussian, with exponential tails in many cases. 
Exponential tails have indeed been observed and can be theoretically ex plained 
by the intermittency of the velocity distribution l|Falgarone et allll994h . 

The linewidth-size relationship found for the ROCs is very flat and does 
not correspo nd to that of Larson 's law described in § 15.11 This has been fully 
discussed bv lOstriker et al.1 l)200l|) . who determined that the least squares fit is 
associated with the superposition of density structure along the line of sight. 
The relation for actual clumps or cores is related to the lower envelope of the 
linewidth-size scatter plot, which has a slope of ~ 0.3 ± 0.1, toward the shallow 
end of the observed values given by equation 1341 

A similar analysis of synthetic CO j (20-19) maps emphasises the differences 
between the isothermal and non-isothermal regimes. Emission maps in this CO 
line are drastically different, as displayed in Fig. The peak intensity in the 
isothermal case is a factor of more than 10 3 lower. This is to be expected, 
as the J = 20 — > 19 transition traces only those regions where the gas has 
T « 1000 K. In the case of MT, hot regions are naturally associated with strong 
shocks, whereas in the (quasi-) isothermal case the temperature is much lower 
and uniform. Therefore, in the isothermal case, the emission in this line simply 
reflects the column density (compare the lower maps of Figs. 1131 and 118(1 . 

For the J = 20 — > 19 emission maps, the values of parameters ()41JI — 1(44(1 
are again given by least-square fits. Emission maps based on the MT simulation 
data yield 



b = (4.7 ± 0.4) x 10~ 2 , b' = (2.7 ± 0.2) x 10~ 2 , 
c = -(5.01 ± 0.04) x 10 _1 , d = 1.950 ± 0.002. 
and for the quasi-isothermal case, 

b = (6.2 ± 0.7) x 10~ 2 , b' = (2.7 ± 0.3) x 10~ 2 , 
c = -(5.71 ± 0.06) x 10" 1 , d = 1.936 ± 0.007. 



(48) 



(49) 



. It is remarkable that statistical properties of the maps of CO(20-19) as sum- 
marised by equations (|48|l and (|49|l do not reveal any differences between MT 
and IT, although the underlying maps are completely different. They are very 
similar to the parameters of CO (4-3) emission, see Eqs. 1(45(1 - 1(47(1 . Linewidth- 
size relationship is somewhat steeper for the CO (20-19) maps, which might 
be related to the size of the statistical sample (number of ROCs is smaller for 
J = 20 — * 19 transition). 

The J — 20 — + 19 emission line profiles derived from molecular data are 
somewhat different from the line profiles derived from the quasi-isothermal data, 
with the emission lines often more prominently peaked (though far weaker) in 
the quasi-isothermal case. 

The remaining parameters are qualitatively similar to the parameters derived 
for the J = 4 — > 3 line emission maps. The virial parameter a is large (a > 
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10 10 ), confirming that the kinetic energy is the dominant form of energy in this 
turbulence regime, and the 'cores' on the maps (ROCs) are false cores which 
cannot be gravitationally confined. The fact that c < indicates that gravity 
is relatively more important for regions with small mass (i.e. of smaller size). 

The scale dependence of ROC mass and the virial parameter a have the 
following interpretation. In the limit of constant column density, the mass of 
a ROC of size A is M oc p co \\ 2 , where p co \ is the column density, and the 
resulting scale dependence is M oc A 2 . These are close to the values we deduce 
fro m all emission maps, M oc A d , where d = 1.96 ± 0.03, as concluded also 
bv iBallesteros-Paredes and Mac Low! l|2002j) . Scaling of the virial parameter a 
follows from the result for the mass scaling. Given that a is independent of A, 
the scaling dependence of the parameter a approaches oc M~ - 5 (in the limit of 
constant column densities). The numerical value is close to this limit: we find 
a oc M c , where c = -0.5 ± 0.2. 

5.4 Discussion 

For molecular line emission that traces the bulk of the gas, isothermal simula- 
tions can be very accurate in mapping molecular emission. Properties of the 
emission maps and shapes of emission lines reproduce the corresponding values 
deduced from the detailed molecular simulations. An attempt to map shocked 
regions with a temperature-sensitive chemical tracer fails in the case of isother- 
mal turbulence, as is naturally anticipated. Results of isothermal turbulence 
simulations should only be used for what they have been made — tracing gas 
close to the average temperature. 

The predictions of values of scaling exponents from analysis of emission maps 
with different CO rotational number J are approximately equal. 

The algorithm for construction of synthetic emission maps used here is very 
simple. It functions only for the case of optically thin emission and may serve 
as a starting point for a deeper analysis. Self-consistent three-dimensional ra- 
diative transfer algo rithms are ess e ntial for obtaining reliable predictions as, 
fo r examp l e, use d bv lPadoan et alJ l)l998[) based on the radiati ve transfer code 
of l.Iuvelal Jl997h . The large velo city gradient approximat i on jSobolevt Il957|) 
for radiative transfer was used by lOssenkopf and Mac Low! l|2002fl to construct 
synthetic emission maps. 

6 CONCLUSIONS 

We have examined how molecular dynamics influences supersonic turbulence 
and, vice versa, how turbulence influences molecular content. Numerical sim- 
ulations of cloud dynamics have usually employed an isothermal equation of 
state as an approximation for gas inside molecular clouds. Hence, we have com- 
pared our results (MT) to that of the equivalent isothermal simulation (IT). We 
analysed decaying turbulence, beginning from a molecular state in which strong 
shocks develop and destroy 87% of the molecules and running until 80% of the 
molecules have reformed. Our main results are as follows. 

1. There is no significant difference between the evolution of power spectra 
for MT and IT. The spectra are not power laws but involve a faster decay 
in the energy associated with the higher wavenumber regime. 
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Figure 18: Maps of 12 CO J = 20 — ► 19 emission constructed from the MT data 
(upper panel). The same emission map constructed from the quasi-isothermal 
turbulence data (lower panel) as described in § 15.21 Note the huge difference in 
intensity as given by the colour bars. 
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2. Despite strongly supersonic root-mean-square velocities, the ratio of en- 
ergy in compressional modes to energy in solenoidal modes at low wavenum- 
bers drops from the initial value of 1:2 to ~1:5 almost immediately. Little 
further decrease occurs at the low wavenumber end while the flow remains 
supersonic. 

3. Compression waves steepen much faster in the simulation of molecular 
turbulence leading to a rapid increase in high wavenumber compressional 
energy. This appears to be due to the slower dissipation of thermal energy 
during the dissociative phase than in the equivalent isothermal phase. In 
contrast, the high wavenumber compressional energy in IT later overtakes 
that of the MT. 

4. The effective poly tropic index in MT varies considerably. In the very 
early stages, there are effectively two isothermal phases, corresponding to 
atomic gas at ~ 8000 K and molecular gas at ~ 40 K. The index is then 
reduced to below and above unity for extended periods (see Fig. 3$ before 
approaching unity at late times. The sub- isothermal period is related to 
the presence of strong compression and molecule reformation behind fast 
shocks. 

5. The density PDFs in the MT case are found to be close to log-normal. 
The strongest deviations occur early on when strong shocks distort the 
statistics. The deviations are smaller in the MT case than in IT, probably 
due to the higher effective polytropic index. 

6. The structure of the velocity fields was studied with the help of the di- 
vergence and modulus of curl. Vortex isosurfaces take the form of sheets 
closely related in space and size to the shock surfaces. 

7. The PDFs of vorticity, as well as velocity gradients, take the form of 
exponentials. We constructed PDFs of compensated distributions where 
the change of the mean values of the distributions is accounted for in the 
exponentials. With constant parameters, the functions provide a useful 
tool for diagnosing decaying turbulence. Velocity PDFs maintain Gaussian 
distributions. 

8. Simulated maps of CO rotational emission lines are constructed. We find 
that isothermal simulations provide excellent agreement with the molec- 
ular predictions for tracers of the cool bulk of the gas. This is confirmed 
through the 'Region of Contrast' method to characterise the properties 
of the emission maps. In Paper I, we showed that as turbulence decays, 
molecules do not remain within swept-up shells but are distributed very 
homogeneously, and the molecular fraction distribution has a small vari- 
ation. This results in very similar looking mass and molecular mass dis- 
tributions during the reformation phase. On the basis of this result it is 
argued that isothermal simulations can be used to map molecules. This 
hypothesis is supported by the synthetic maps. 

9. On the other hand, large differences are found in the emission distribution 
of high-J CO lines. These lines are sensitive to the temperature. 
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10. In ast rophysical turbulenc e, the magnetic field influences turbulent cas- 
cades ijVestuto et allEoO^ . For example, the column density power spec- 
trum was found to be significantly shallower as the field strength is raised 
llPadoan et all Eool . A strong magnetic field also enhances the shock 
number transverse t o the field direction at the expense of parallel shocks 
ijSmith et all 12000^1 . 

The magnetic field would also weaken shocks, allowing molecules to survive 
the passage of stronger shocks. This will tend to reduce the differences 
with respect to the IT case caused by the chemical and cooling processes. 

The main conclusion of this work is that isothermal simulations adequately 
model molecular turbulence and can be used to predict many properties of the 
molecular emission; given the supersonic dynamics, the chemical reactions in 
a turbulent gas can be significantly accelerated due to strong compression and 
advection. 

This work provides a basis for future research in a number of directions. It 
is possible that the behaviour of molecules can be modelled as a passive scalar, 
and pseudo-temperature distributions can be constructed from appropriately 
gauged divergence fields. 
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